Cold dilute neutron matter on the lattice I: Lattice virial 
coefficients and large scattering lengths 

Dean Lee and Thomas Schafer 
Department of Physics, North Carolina State University, Raleigh, NC 27695 

Abstract 

We study cold dilute neutron matter on the lattice using an effective field theory. We work in 
the unitary limit in which the scattering length is much larger than the interparticle spacing. In 
this paper we focus on the equation of state at temperatures above the Fermi temperature and 
compare lattice simulations to the virial expansion on the lattice and in the continuum. We find 
that in the unitary limit lattice discretization errors in the second virial coefficient are significantly 
enhanced. As a consequence the equation of state does not show the universal scaling behavior 
expected in the unitary limit. We suggest that scaling can be improved by tuning the second virial 
coefficient rather than the scattering length. 



I. INTRODUCTION 



Cold dilute neutron matter is an intriguing physical system. It is relevant to the physics 
of the inner crust of neutron stars 1]. It is also close to an interesting universal limit which 
is the result of a large separation of scales. The neutron scattering length is a nn ~ —18 
fm, which implies that the dimensionless parameter kp\a nn \ 3> 1 for densities p > 10 _4 pjv. 
Here, kp = (37r 2 p) 1,/3 is the Fermi momentum and p^ ~ 0.16 fm -3 is the saturation density 
of nuclear matter. The effective range, on the other hand, is r nn ~ 2.8 fm. So if the 
density is very small, p < O.lpN, then kp\r nn \ is a small parameter and neutron matter 
is close to the limit in which kp \a nn \ — > oo and kp\r nn \ — > 0. This is known as the 
unitary limit, where the vacuum scattering amplitude in the s-wave channel has a zero- 
energy resonance, and the cross-section saturates the unitarity bound. In this limit there 
is no expansion parameter and the calculation of the equation of state and of transport 
properties is a difficult non-perturbative problem. It is now possible to create systems 
of fermions in the unitary limit in the laboratory by trapping fermionic atoms and tuning 
the scattering length using a Feshbach resonance . This technique will provide 

experimental measurements of universal parameters, but it is clearly desirable to also develop 
a computational approach. There are several computational studies of neutron matter at 
zero temperature using potential models and Green's function Monte Carlo [?], 0, H liot . 
Recently, there have also been simulations on the lattice using effective field theory ll|, [l2j, 
HUG. The main advantage of the effective field theory/lattice approach for the dilute 
neutron matter problem is that it is not restricted to zero temperature. 

We have described our method and some initial results in [l^ . This is the first in a 
sequence of two papers in which we perform a careful study of the parameter dependence 
and scaling behavior of the results. We are particularly interested in verifying that the 
lattice results satisfy the scaling relations that are expected to hold in the unitary limit. 
The main focus of this first paper is to understand the equation of state at low density 
and high temperature and compare the results to the virial expansion. We find that if the 
scattering length is large lattice discretization errors can be as large as 100% or more. At 
fixed lattice spacing when the temperature decreases we find that the error first increases 
before eventually decreasing. In order to counter this effect, we propose that the interaction 
coefficient be tuned to give the correct second order virial coefficient, b2{T), at the chosen 
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FIG. 1: Diagrams contributing to neutron-neutron scattering. 



simulation temperature. We describe in detail how this is done. In the sequel to this paper 
we apply this technique to simulations of cold dilute neutron matter in the unitary limit 
and study the scaling behavior of the results. 

II. LATTICE ACTION 

We consider the theory defined by the partition function 

Z G = Tr exp \-(3{H - pN)] ~ z e 2 ^ L3 J DsDc'Dc* exp [-S] , (1) 
where the lattice action is given by 

S = [ e_/ "X(™K(™ + 6) - e v /3 c^ S (n)+^( 1 _ 6/i)c*(nK(n) 

n,i 

-hJ2 [cM«H)+c;(«K(«4)] +\j2 s2 ^- ( 2 ) 



n.L.i 



Here, n labels the sites of a 3 + 1 dimensional lattice, l s (s = 1, 2, 3) is a spatial unit vector, 
and i labels the two spin components of the neutron, f and j. The spatial lattice spacing 
is a and a t = a t /a is the ratio of the temporal to the spatial lattice spacing. The spatial 
volume of the lattice is L 3 and the temporal length is (3 = 1/T. Dimensional quantities like 
the chemical potential [i and the nucleon mass are given in units of the lattice spacing 
a. We have also defined h = a t /(2mjv)- The Grassmann fields are denoted by Q(n) and 
s(n) is a Hubbard-Stratonovich field. 

The interaction coefficient C must be determined for given lattice spacings a and at- In 
jl^ this was done by adjusting C to reproduce the correct neutron scattering length at zero 
temperature. This requires summing the two-particle scattering diagrams shown in Fig. d 
The pole in the scattering amplitude is then compared with Luscher's formula for energy 
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levels in a finite periodic box |16l 1 171]. 

_ 47ra scat t 
m^L 3 

where c x = -2.837297, c 2 = 6.375183. 
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III. RESULTS FOR LARGE SCATTERING LENGTHS 



We present lattice simulation results for the energy per particle as a function of density 
for several different scattering lengths. We are interested in cold dilute neutron matter 
where all length scales are much larger than the lattice spacing. We use a spatial lattice 
spacing of a = (50 MeV) -1 and temporal lattice spacing of a t = (24 MeV)- 1 . We will 
consider scattering lengths |a sca tt| > 4.675 fm. The spatial lattice spacing was chosen so 
that it is smaller than the smallest scattering length. The temporal lattice spacing was 
chosen sufficiently small so that the results are very close to the a t — > limit. We have 
computed the energy per neutron for five different scattering lengths. The corresponding 
operator coefficients are shown in Table 1. 

Table 1: Operator coefficients and scattering lengths 



C (10- 4 MeV- 2 ) 


-1.028 


-1.108 


-1.153 


-1.257 


-1.318 


(10" 5 MeV" 2 ) 


-2.259 


-2.239 


-2.218 


-2.141 


-2.078 


a scatt (fm) 


-4.675 


-9.35 


-18.70 


+18.70 


+9.35 



We have also shown the derivative jR- which is needed to compute derivatives with respect 

t0 I use t be ^ Moabe Co Wm Q to „ Hubbard-Stratonovich Md 
configurations as described in [l^. We use diagonal preconditioning before each conjugate 
gradient solve. If K is the single-spin neutron matrix for a given Hubbard-Stratonovich 
field configuration, then we must solve the linear equation 

K ] Kv 



b. 

Rather than solving this directly, we make use of the diagonal matrix 

D = diag [K*K] 

and solve instead 

[D^K^KD- 1 ] Dv = D- X b. 
4 



(4) 



(5) 



(6) 



Fluctuations associated with the Hubbard-Stratonovich field occur on the diagonal of K, 
and therefore the matrix D^ 1 K> K D^ 1 typically has a smaller condition number than K^K. 
We further improve the performance by adding a small positive constant e to produce the 
modified equation 

[D^K^KD- 1 + e] Dv = D~ x b. (7) 

We tune e so that the equilibration time for the hybrid Monte Carlo algorithm is minimized 
while keeping the induced systematic error smaller than the stochastic error. In practice 
we take e to be 10~ 4 or smaller. 

Roughly 10 4 HMC trajectories were run, split across 4 processors running completely 
independent trajectories. Averages and errors were computed by comparing the results of 
each processor. The finite volume error was tested by going to larger volumes. The final 
lattice sizes were chosen so that the finite volume error was less than one percent. For the 
data at T = 4 MeV we used a lattice of size 4 3 x 6. For T = 3 MeV we used 5 3 x 8, and for 
T = 2 MeV we used 5 3 x 12. As an example of the finite volume dependence, we show in 
Table 2 the density and energy per particle from simulations with lattice volumes 4 3 , 5 3 , 6 3 
at a scatt = —18.70 fm, T = 4 MeV, and \x = 0. For comparison we also show the volume 
dependence of the bubble chain calculations described in [13j and Sect. IVl 



Table 2: L dependence for a scatt = —18.70 fm, T = 4 MeV, /i = 



L 


p bubblc ( 10 -3 fm -3) 


^bubble (MeV) 


p sim ( 1Q -3 fm -3) 


f^(MeV) 


4 


7.688 


3.527 


8.70(3) 


3.41(3) 


5 


7.695 


3.521 


8.78(10) 


3.41(3) 


6 


7.695 


3.521 


8.81(10) 


3.39(3) 



These results suggest that the finite volume effects are smaller than the statistical errors. 
If we take the volume dependence from the bubble chain calculations as a guide, then the 
finite volume errors are well below one percent. 

Fig. 121 shows the lattice simulation results for the energy per neutron at T = 4 MeV. Fig. 
El shows the energy per neutron at T = 3 MeV, and Fig. 0] shows the energy per neutron 
at T = 2 MeV. In each of the plots we have taken a range of densities from zero to a 
quarter-filled lattice. With a spatial lattice spacing of (50 MeV) -1 , the quarter- filled lattice 
corresponds with a density of 0.0081 fm -3 . Beyond this one might find significant lattice 
artifacts. We observe that the energy per particle depends quite strongly on the scattering 
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FIG. 2: Energy per neutron versus density at T = 4 MeV for various scattering lengths. 
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FIG. 3: Energy per neutron versus density at T = 3 MeV for various scattering lengths. 



length, both at large and at small density. While the dependence of E/A on the scattering 
length for a degenerate Fermi gas is a complicated, non-perturbative problem it is possible 
to compare our results to theoretical predictions in the opposite limit of low density and 
high temperature. 
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FIG. 4: Energy per neutron versus density at T = 2 MeV for various scattering lengths. 



IV. VIRIAL EXPANSION 

The virial expansion arranges multi-particle interactions as a power series in fugacity, 

z = e^. (8) 

For example the logarithm of the partition function per unit volume can be written as 

i In Z G = PP = A [z + b 2 (T)z 2 + b 3 (T)z 3 ■■■]. (9) 

The second order virial coefficient b 2 (T) is determined entirely by two-particle interactions, 
while the third virial coefficient depends on three-body interactions, and so on. We can use 
the virial expansion to compute thermodynamic observables in the high-temperature/low- 



density limit. The virial expansion is reliable if the thermal wavelength \ T = a/ (2ty) / (toatT) 
is smaller than the interparticle spacing, At < p" 1 / 3 . The neutron density can be computed 
in terms of the derivative of In Z G with respect to the chemical potential, 

" = ^wi lnZa - (10) 

To second order in the virial expansion we find 

P = ^[z + 2b 2 (T)z 2 + •••]. (11) 
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As a result of discretization error we find on the lattice a more general power series in 
fugacity, 

p = ^b 1 (T)[z + 2b 2 (T)z 2 + -..]. (12) 

While b\(T) is no longer guaranteed to equal 1, we should find that b\(T) — ► 1 in the 
continuum limit. We will use (|T2*|) as the definition for the virial coefficients on the lattice. 

In the unitary limit, where the effective range is zero and scattering length is infinite, one 
finds 3 

6 2 (T) = 3 ■ 2 _ 3 « 0.530. (13) 
We can compare this with the result for a free Fermi gas, where 

b 2 {T) = -2"3 « -0.177. (14) 

For zero range but finite scattering length the second virial coefficient becomes 

f ^[l-erf(N)]-^ forx<0, 

6 2 (T) = < , Bfl | (15) 
|^ v/2e fe BT _ e_ [x _ erf(x)] _ _i_ for x > o, 

where erf is the error function, E B is the two-particle bound state energy for positive scat- 
tering length, and 

x = -=t^ . (16) 

V 2vra sca tt 

As the effective range goes to zero we have the relation 

\ E b\ = z4— , (17) 



ma scatt 



and therefore we can write 



^[l-erfdxl)]-^ forx<0, 



b 2 (T) = { 2 (18) 

V2e* 2 - 7j [1 - erf {x)\ - ^ for x > 0. 

See |20j for a calculation of b 2 (T) with realistic interactions, taking into account effective 
range corrections as well as higher partial waves. 



V. BUBBLE CHAIN DIAGRAMS AND THE VIRIAL EXPANSION ON THE 
LATTICE 

The second virial coefficient is determined by two body interactions. This means that 
we do not have to rely on simulations in order to determine b 2 (T) on the lattice, but we can 
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FIG. 5: Bubble chain diagrams contributing to the neutron self-energy. 



extract the second virial coefficient from the lattice regularized bubble chain diagrams. For 
the neutron propagator we have the diagrams shown in Fig. 03 The bubble diagrams give a 
contribution to neutron self-energy of the form 

D iree (p- q) 



:i-^(e--i)^-ET37r 



6/z) 2 (e- 



-Ca t 



(19) 



where 
B 



L 3 L t ^ e-^te-iP'o/ie-^* - 1 + u) p / 2 +k e'^e-^o^e^* - 1 + U- 
We use this to compute the full neutron propagator 



p/2+k 



(20) 



D iul \q) 



The average number of neutrons is then 



(3 o/j, 



l-j:(q)D hcc (q)' 



{m N -n)a t ^ 



(21) 



(22) 



For the logarithm of the grand canonical partition function, InZc, the lowest non-trivial 
order in pA^ is given by the bubble chain diagrams shown in Fig. El These give a contribution 



In Z G - In Z G 



free 



1 -In [1 - (1 -6h) 2 (e~ Ca 



1) B(p + q,fj,)] D ivcc (p)D ivee (q) 



p.q 



B(p + q,fi) 



(23) 



We can now compare these results to the results from simulations and to the virial 
expansion in the continuum. In Fig. Owe plot the density versus fugacity at T = 4 MeV 
for scattering lengths a scat t = ±18.7 fm. We show data for a free Fermi gas on the lattice, 
bubble chain calculations, and full simulation results. We also plot first and second order 



FIG. 6: Bubble chain diagrams contributing to the logarithm of the partition function. 
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FIG. 7: Density versus fugacity at T = 4 MeV for scattering lengths a sca tt = ±18.70 fm. We show- 
comparisons with the first and second order virial expansion. 

virial results using &i(T) = 1.27, which we obtained by fitting the bubble chain data at 
very small fugacity. The results for the free Fermi gas on the lattice agree with the second 
order virial results with 6 2 [T) — —0.177. The bubble chain and full simulation results for 
Oscatt = ±18.70 fm are not very far from the second order curve for bz(T) = 0.530. However 
the results for —18.70 fm and ±18.70 fm both lie above the second order virial curve, rather 
than one on either side. 

In Fig. |H] we plot the density versus fugacity at T = 3 MeV and virial curves with 
bi(T) = 1.28. In Fig. El we plot the density versus fugacity at T = 2 MeV and virial curves 
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FIG. 8: Density versus fugacity at T = 3 MeV for scattering lengths a sca tt = ±18.70 fm. We show 
comparisons with the first and second order virial expansion. 



with bi(T) = 1.22. Again we determined bi(T) by fitting the bubble chain data at very 
small fugacity. The free Fermi gas results on the lattice at T = 3 MeV and 2 MeV match 
the second order virial results at b 2 (T) = —0.177. But we find the same problem for the 
bubble chain and full simulation results. The results for a scat t = —18.70 fm and +18.70 fm 
both lie above the second order virial curve for b 2 (T) = 0.530. Furthermore the deviation 
appears to be worse at colder temperatures. 



VI. LATTICE VIRIAL COEFFICIENTS 



We use the bubble chain expansion to compute b 2 (T) on the lattice for the various coupling 
strengths and temperatures mentioned above. The specific procedure we use is as follows. 
We first compute the free Fermi gas density on the lattice and use the relation 
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p frcc w — 6i (T) z - 2 • 2~2 Z 2 | (24) 

Arp 



to determine b±(T). This can be measured at any small fugacity, and we choose z = e~ 5 ss 
0.0067 or /j, = —5k B T. In Table 2 we show results for &i(T) for a range of temperatures. 
We have taken the lattice volume to be large enough so that the finite volume error is 0.1% 
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FIG. 9: Density versus fugacity at T = 2 MeV for scattering lengths a sca tt = ±18.70 fm. We show 
comparisons with the first and second order virial expansion. 



or less, and the values for L are shown in Table 2. 

Table 2: b x {T) on the lattice 



T (MeV) 


L 


h{T) 


4 


5 


1.273 


3 


6 


1.282 


2 


6 


1.221 


1.5 


7 


1.160 


1 


8 


1.091 


0.667 


9 


1.053 


0.5 


10 


1.038 



As the temperature decreases we see that &i(T) — > 1. This is expected since at fixed lattice 
spacing the system moves closer to the continuum limit as we decrease the temperature. 

With the same chemical potential and lattice volumes, we compute the density in the 
bubble chain approximation and determine b 2 (T) using the relation 

p bubbic w * bl (T) [z + 2- b 2 (T)z 2 } . (25) 
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The results for b%(T) for the various coupling strengths are shown in Table 3. For each pair 
shown the lattice result is on the left and the continuum result is on the right. 

Table 3: Second virial coefficient &2CO for different scattering lengths a sca tt as 
calculated on the lattice (first column) and in the continuum limit (second column). 



a scatt (fm) 


-4.675 


-9.35 


-18.70 


+18.70 


+9.35 


C (lCT 4 MeV~ 2 ) 


-1.028 


-1.108 


-1.153 


-1.257 


-1.318 


T = 4 MeV 


0.66 


0.198 


0.88 


0.322 


1.03 


0.411 


1.45 


0.692 


1.77 


0.917 


T = 3 MeV 


0.76 


0.170 


1.07 


0.299 


1.29 


0.396 


1.97 


0.722 


2.52 


1.004 


T = 2 MeV 


0.80 


0.131 


1.26 


0.263 


1.61 


0.371 


2.85 


0.776 


3.98 


1.176 


T = 1.5 MeV 


0.71 


0.103 


1.23 


0.237 


1.68 


0.352 


3.42 


0.825 


5.22 


1.350 


T = 1 MeV 


0.46 


0.066 


1.09 


0.198 


1.44 


0.322 


3.79 


0.917 


6.85 


1.720 


T = 0.667 MeV 


0.23 


0.031 


0.60 


0.159 


1.01 


0.289 


3.58 


1.047 


8.01 


2.368 


T = 0.5 MeV 


0.12 


0.008 


0.40 


0.131 


0.74 


0.263 


3.29 


1.176 


8.82 


3.167 



We see that the lattice virial coefficients are larger than the continuum values. Also as we 
decrease the temperature at fixed scattering length, the deviation in 6 2 {T) first increases 
before eventually decreasing. 



VII. DISCUSSION 

In Fig. E3we plot &2CO versus inverse scattering length for a range of temperatures. The 
smooth curves are the continuum results given in ()18j) . and the points with interpolating 
guide lines are the lattice results from Table 3. For the continuum curves we see that the 
value of 62 at infinite scattering length remains fixed at 3 • 2^2 « 0.530. However as the 
temperature decreases, the slope as a function of inverse scattering length steepens. The 
functional dependence of 6 2 on a scat t and T is roughly 

b 2 ~ e x2 = e a ™* m " T . (26) 

We observe that the lattice virial curves are actually not very different from the continuum 
virial curves. However they are shifted to the left slightly as a function of the inverse 
scattering length. The large slope together with this leftward shift is responsible for the 
large lattice virial coefficients. This is seen even more clearly in Figs. ^2 and El where 
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FIG. 10: Plot of second virial coefficient as function of inverse scattering length for various tem- 
peratures. The smooth curves are the continuum results, and the points with interpolating lines 
are the lattice results. 

we show the continuum and lattice virial results at temperatures T = 1 MeV and 0.5 MeV 
along with shifted lattice results so that 62 equals 0.530 at infinite scattering length. We 
see that the shifted curves appears to match the continuum results relatively well. The 
required shift as a function of temperature appears to decrease with temperature, though 
the exact dependence requires further study. 

Physically the large error in the second virial coefficient is directly related to singularities 
in the sum over particle-particle bubble diagrams due to either true bound states for positive 
scattering length or zero-energy resonances in the unitary limit. Although the singularities 
are completely integrable at nonzero temperature, on the lattice the momentum integrals 
are finite sums which can be very sensitive to small perturbations. The obvious way to 
address this problem, aside from working at a much smaller lattice spacing, is to use an 
improved lattice action, one that includes more than simple nearest-neighbor hopping terms 
for the kinetic energy as well as higher order corrections to the interaction. We are currently 
investigating improved actions, but this is not an easy task and not all improved operators 
maintain the positivity of the euclidean action. 

We propose that even with an unimproved action one can improve the scaling behavior 
of the data by tuning the coefficient C to give the correct physical value of the second virial 
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FIG. 11: Plot of continuum and lattice virial results for T = 1 MeV. We also show the lattice 
results shifted horizontally so that b 2 equals 0.530 at infinite scattering length. 




T= 0.5 MeV 
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FIG. 12: Plot of continuum and lattice virial results for T = 0.5 MeV. We also show the lattice 
results shifted horizontally so that b 2 equals 0.530 at infinite scattering length. 

coeffcient b 2 (T) for each desired simulation temperature. Since b 2 (T) on the lattice is easily 
calculated this is not much of an added burden. This procedure coincides with the standard 
approach of fixing the zero temperature scattering length in the limit that the lattice spacing 
goes to zero or as T — > at fixed lattice spacing. Fixing b 2 (T) will obviously improve the 



15 



scaling behavior at low density. Whether it will also improve the data in the degenerate 
regime is not a priori clear and this question will be the focus of our companion paper. 
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